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, O-f - Abstract 

| We study synchronization of locally coupled noisy phase oscillators that move diffusively in 

> 

' a one-dimensional ring. Together with the disordered and the globally synchronized states, the 

IT> ' 

CN ■ system also exhibits wave-like states displaying local order. We use a statistical description valid 

. for a large number of oscillators to show that for any finite system there is a critical mobility above 

O 

which all wave-like solutions become unstable. Through Langevin simulations, we show that the 
transition to global synchronization is mediated by a shift in the relative size of attractor basins 
^ associated to wave-like states. Mobility disrupts these states and paves the way for the system to 

h : 

attain global synchronization. 
PACS numbers: 05.45.Xt,87.18.Gh 
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Synchronization of oscillators is a widespread phenomenon in nature jl-3]. In biology, 
synchronization can occur at scales that range from groups of single cells to ensembles of 
complex organisms [4]. When oscillators hold fixed positions in space and the interaction that 
drives synchronization is short ranged, spatial and temporal patterns can self-organize. Such 
is the case in cardiac tissue, where cells generate spiral patterns that shape the heartbeats j^J. 
Also in central pattern generators, the oscillating neural network self-organizes to produce 



coordinated movements of the body 



6]. 



A different situation arises when the oscillators are not fixed in space but are able to 
move around. The problem of synchronization of moving oscillators has many applications 
in the domain of chemistry 3], biology [8], and technology [9]. Small porous particles 
loaded with the catalyst of the Belousov-Zhabotinsky reaction behave as individual chemical 
oscillators, undergoing a density-dependent synchronization transition as the stirring rate is 
increased The same particles support wave propagation in the form of dynamic target 
and spiral patterns when the particles are not moving [h]]. This phenomenon illustrates a 
wider scenario: mobility and mixing remove local defects and patterns, enabling global order. 
This effect has far reaching consequences in finite systems. For example, in ecosystems of 
competing populations with cyclic interactions, biodiversity can be sustained if dispersal 
is local, but it is lost when dispersal occurs over large length scales ll||. The dynamics 
of such cyclic competition was described by a complex Ginzburg-Landau equation near a 
Hopf bifurcation, displaying complex oscillatory patterns indicative of biodiversity for low 
mobility, while in the case of high mobility diversity is wiped out [12I H| • 

In this paper we study the effects of mobility — spatial diffusion — on the macroscopic 
collective dynamics of locally coupled, moving phase oscillators subjected to noise, in a one- 
dimensional ring. When oscillators are fixed in space, these systems can exhibit a series 
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16fl . Such states have been called m- twist 



of steady states where local order is present 

solutions [lS |, Fig. [TJ Here we show that mobility can destabilize all m- twist solutions, 
enhancing the stability of the synchronized solution. We find that in finite systems there is 
a critical mobility above which either the synchronized or the disordered state is stable. 
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I. DIFFUSING PHASE OSCILLATORS 



We consider an ensemble of N identical phase oscillators that diffuse on a ring of perime- 
ter L. Oscillators are coupled to other oscillators in their local neighborhood, within an 
interaction range r. The dynamics of phase and position is described by 
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where i — 1, . . . , N is the oscillator label, 8i(t) and Xi(t) are the phase and position of the 
i-th oscillator at time t, oj is the autonomous frequency, and 7 is the coupling strength 
— whose inverse characterizes the typical relaxation time of the interaction. Each oscillator 
interacts with its neighbors in the range r through the coupling function in brackets, which 
defines an attractive interaction towards the local average of the phase. In steady state the 
spatial density is uniform and the number of neighbors is on average constant, n — = 
Nr/L. With this definition of the coupling the thermodynamic limit is well defined, and 
the system reduces to the noisy Kuramoto model for r = L/2 [l7]. The fluctuation terms 
£o t i and ^ X)i represent two uncorrelated Gaussian noises such that (Ze,i(f)) = (6m W) = 0' 
and (£0,1 (t)^,.,-^')) = ((,x,i(t)£, X j(t')) = 5ij5(t — t'). The strength of angular fluctuations is 
determined by the angular diffusion coefficient C, while the spatial diffusion coefficient D 
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FIG. 1: The system, Eqs. {U-©, exhibits different kinds of states: (a) disorder, (b) partial 
synchronization, and twisted states: e.g. (c) 1-twist and (d) 2-twist. Each dot (x, 9) represents 
an oscillator. Parameters in arbitrary units: N = 1000, L = 2ir, r = L/100, oj = 0, 7 = 0.1 and 
V2D = 0.01. In (a) V2C = 0.35 and in (b)-(d) V2C = 0.10. States (b)-(d) coexist. The snapshots 
were taken after 15000 time units, starting from random initial conditions. 
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FIG. 2: Twisted states affect the ensemble average of the global order parameter, (Z). (a) (Z) vs. 
scaled angular diffusion C/7 for D/7 = 0.03 (green squares) and D/7 = 0.39 (blue dots). For small 
D/7, (Z) does not reach one, even for C = 0. Vertical dotted lines correspond to values of C/7 
in (b). (b) (Z) vs. scaled mobility -D/7, for C/7 = 0.11 (orange triangles) and C/7 = 0.25 (red 
diamonds). Notice that as D is reduced, (Z) decays. Vertical dotted lines correspond to values 
of C/7 in (a), (c) The histogram of Z splits into two different peaks at low and high values of Z 
for low values of -D/7, but displays only one peak at high Z for high D/7. Other parameters are 
N = 1000, L = 2vr, r/L = 0.01, u = and 7 = 0.1. 

determines the mobility of oscillators. The oscillators are point-like particles that can diffuse 
freely, i. e. their movement is not affected by the presence of other oscillators. Both the phase 
6i(t) and position Xi(t) are periodic variables such that < 9i(t) < 2n and < Xi(t) < L. 



II. GLOBAL ORDER PARAMETER 

The system described by Eqs. ([I])-© displays a range of states illustrated in Fig. [TJ 
We can characterize global order by the order parameter Z, defined as the absolute value 
of the population average of the complex unit vectors of the oscillator phases Z(t) = 
N^ 1 ] e* 61 ^! [18]. Mobility introduces a dramatic change in the ensemble average of 
the order parameter (Z), where (...) denotes an average over different initial conditions 
and realizations of the noise, Fig. [2j When mobility is large enough, the system displays 
global synchronization as C — > 0, blue dots in Fig. [2][a). However, when mobility is reduced, 
global order is compromised, as indicated by the green squares in Fig. Ufa). Increasing mo- 
bility D results in an increase of the ensemble average (Z), Fig. W(Jd)- The cause for this 
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behavior can be traced back to the existence of twisted states, which display local order 
but have themselves a vanishing global order parameter, Fig. QJc,d). Due to the presence 
of these states, (Z) results from the average of a bimodal distribution P(Z), with a peak 
related to the twisted states and the other to global synchronization. Above a critical value 
of the mobility, P(Z) becomes unimodal, Fig. Thus, although a global parameter is 

not suited to capture the complexity of a system with local interactions, its statistics reflect 
the existence of twisted states. 



III. STATISTICAL DESCRIPTION 



The role of twisted states can be studied using a statistical description that is valid when 
the number of oscillators is large. Given that the oscillators have identical autonomous 
frequencies u, it is convenient to make the transformation 9 — > 9 — cut to a rotating reference 
frame. We coarse grain the microscopic model and describe the system in terms of p(x, 9, t), 
the density of oscillators at position x with phase 9, which obeys the Fokker-Planck equation 

d t p(x,9,t) = Dd xx p(x,9,t) + Cd ee p(x,9,t) (3) 

fL r 2ir 



nix) 



nL/ PZ7T 

/ dx' de'g{x-x')sm{e -e , )p{x , ,e , ,t)p{x,e,t) 

.Jo Jo 



where g(x—x') is a kernel accounting for the range and relative strength of local interactions, 
while 

n{x)= dx' d9'g(x-x')p(x',9',t) 



JO JO 

denotes the effective number of oscillators in this range. In this paper we choose g(x—x') = 1 
for | a; — x'\ < r and g(x — x') — otherwise, as in Eq. (TjQ). The derivation of Eq. ([3]) relies 
on the assumption that p 2 (x,9,t; x' ,9' ,t) = p(x,9,t)p(x' ,9' \t) [3]. 

Since the movement of the oscillators is purely diffusive, see Eq. (T5]), the spatial density 
of oscillators is uniform in steady state, J Q 27r d9p(x,9,t) = N/L = po, and n(x) = 2rpo- For 
small N, fluctuations in the spatial density can induce the formation of gaps in which the 
nearest oscillator is beyond the range of interaction. In this paper we consider large densities 
such that the lifetime of these gaps is much shorter than other typical time-scales. 



5 



A. Local order parameter 



The statistical description fl3]) can be cast in a more transparent form introducing a local 
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mean field. Local order can be characterized by a local order parameter 

R(x, t)j*W = tdx' fde' 9 ^'^ e i0 'p(x', 9', t) , (4) 
Jo Jo n(x) 

where R(x, t) is a measure of local order and ip(x, t) is the local average of the phase. Eq. (J3]) 
can be expressed in terms of this local order parameter as 



d t p{x, 9, t) = Dd xx p(x, 9, t) + Cd ee p(x, 9, t) 

+ 7i2(ar, t)d e (sin (9 - ip(x, t)) p(x, 9, t)) 



(5) 



reflecting the fact that if)(x, t) acts as a local mean field and R(x, t) is a local modulation to 
the coupling strength. 



IV. TRANSITION FROM DISORDER TO LOCAL ORDER 

Eq. 03]) has a trivial steady state p(x,9,t) = po/2n = p d which corresponds to the 
disordered state of the system. We study the stability of pd bv inserting p(x,9,t) = pa + 
ef(x,t) cos (£9) in Eq. ([3]) and keeping terms of order 0(e) [20]. Linear stability analysis 
reveals that the disordered solution pd becomes unstable for £ = 1 when C < C*, with 

C* = 7 /2 . (6) 

This threshold is independent of po and D, and determines the value of C below which 
local order sets in. The critical C* given by Eq. ([6]) coincides with the existence [l^] and 
stability [2o| threshold displayed by globally coupled noisy oscillators. 



V. LOCAL ORDER SOLUTIONS 



Once local order has set in, the system also supports twisted solutions. We specifically 
look for steady state solutions to Eq. (jSJ) of the form p s (x,9) = f{9 — ip(x)). Such wave- 
like solutions describe densities in which the angular distribution has the same shape, but 
is centered at position dependent phases ip(x). Setting d t p = we obtain an ordinary 
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differential equation for / 

D e + D x (^'(x)) 2 ) f\<p) + 1 R (sin(^) - D x r(x)) f{<p) = C(x) , (7) 

that we can solve together with periodic boundary conditions in phase and space to determine 
the arbitrary constant C(x) and phases ip(x). Periodicity of the phase is consistent with 
solutions that fulfill ip"(x) = and C(x) = 0, and periodicity of space sets the wave numbers 
k = 2iTm/L with m integer. We obtain the m- twist steady state solutions 



p s (x, 6) = Afexp 



'eff 
2tt 



cos(6 ) — kx) 



where J\f is a normalization constant such that J Q n d6 J Q dxp s (x, 
We have introduced the effective diffusion coefficient 

V eff (k) = C + k 2 D , 



(8) 

iV = Lp , see Fig.^a). 



(9) 



which is a combination of angular diffusion C, and mobility D scaled by the square of the 
wave number k. Effective diffusion competes with the local coupling and controls the 
width of the angular distribution, which has a mean (ip) = kx and variance 

a 2 = 1 - h^R/V^/I^R/Veff) = 1 - R/wac(kr), 



where sinc(x) = sin(x)/x and 



In(z) 



2- 



d9 cos™ 9 exp (z cos ( 



is the modified n-order Bessel function of the first kind. The functional form o 
state angular distribution Eq. (JHJ) is known as the circular normal distribution 



re steady 



2l|. 



A. Effective diffusion controls the local order parameter 

The m-twist solutions described by Eq. (jBJ) exist if and only if the self-consistency condi- 
tion obtained by inserting Eq. (JHJ) into Eq. (TJJ is fulfilled 

flm = si nc(M ^^ , (10) 

see Fig. [3]^b,c). A trivial solution to Eq. ( ITOj) is R m = 0. Apart from this, an expansion of 
Eq. (fTU|) for R m 1 reveals that non-vanishing solutions 

R m « V8V cl - 2 (V c - V eff ) 1/2 (11) 
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FIG. 3: Twisted solutions and local order, (a) Angular distribution of steady state twisted solu- 
tions, Eq. ([8j) . (b) The local order parameter R m decreases with increasing angular fluctuations 
for all solutions, Eq. (|10p . (c) The local order parameter of twisted solutions also decreases with 
mobility, but is not affected for global order, Eq. (|1U|) . The range of interaction is r/L = 0.01. 

exist for T> e g < (7/2) sinc(fcr) = T> c . Again, we see that D e g controls the growth of the local 
order parameter R m characterizing the emergence of twisted solutions, through changes in 
phase fluctuations C and mobility D, Fig. E](b,c). 



B. Existence of twisted solutions 



We can unfold the effects of spatial and angular fluctuations by writing V e ff in terms of 
its components D and C . Setting R m = in Eq. (TTTT) we get 

C + {2TxmjVf D = (7/2) sine (27rmr /L) , (12) 

where we have expressed k in terms of L to stress system size dependence. For each value 
of m and range of interaction r, the surface in parameter space defined by Eq. (|T2l) encloses 
the region where the m-twist solution exists. There are two ways in which fluctuations can 
destroy twisted solutions, by increasing either angular fluctuations or mobility. Mobility gets 
amplified for higher order modes as m 2 , causing twisted solutions to disappear sequentially 
as mobility increases, Fig. HI 

For m = 0, Eq. (TTSl) reduces to C = C* = 7/2 and corresponds to global synchronization. 
Existence of global synchronization in steady state is not affected by mobility, as indicated 
by the dotted red line in Fig. H(a,c). However, existence of twisted solutions in finite systems 
is controlled by angular diffusion C, mobility D, and range of interaction r as indicated by 
Eq. P]), Fig.H 
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FIG. 4: Mobility controls existence of m-twist solutions. Three representative cuts of the phase 
diagram are shown for m = 0,1, 2,4. (a) D = 0, (b) C/7 = 0.1 and (c) r/L = 0.01. The dotted red 
line indicates the onset of local order, which is independent of D. Below the solid green line, the 
dashed blue line and the dot-dashed yellow line, the 1-twist, 2- twist, and 4- twist solution exists, 
respectively, see Eq. (|12j) . An increasing number of co-existing twisted solutions can be found with 
decreasing mobility D and decreasing coupling range r. 



As L — > 00, the critical value of the spatial diffusion coefficient diverges for all m. There- 
fore, in the infinite system size limit all twisted solutions coexist with global synchronization 
for any finite D. These result is in agreement with 22], and indicates that identical noisy 
phase oscillators cannot exhibit a global synchronized state in ID in this limit. 



C. Stability of twisted solutions and states 

While existence and stability thresholds coincide for the global order solution, Eqs. (JH]) 
and (fl2l) . this is not the case for m-twist solutions in finite systems, Fig. [5j We address the 
stability of the 1-twist solution by performing a numerical study of Eq. ([3]), using a finite 
difference scheme. To estimate the stability boundary, we continue a stable twist solution 
until it becomes unstable against small perturbations. We find that the instability is of 
modulational type. The m-twist solutions become stable only after the corresponding local 
order parameter R m becomes larger than a certain value, i. e. twisted solutions become stable 
with a finite amplitude, Fig. |5j For vanishing spatial and angular diffusion C = D = 0, 
we encounter the system studied by Wiley et al. see purple open circle in Fig. E|a). 
The numerical solution seems to approach this point as C/7 — > 0, but the numerics become 
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(a) D = (b) C/Y = 0.1 (c) r/L = 0.01 



disorder 
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FIG. 5: The domains of existence and stability of m-twist solutions do not coincide. The 1-twist 
solution exists below the solid green line, Eq. (fT2|) . but it becomes stable below the dashed purple 
line, determined numerically. The purple dots show the stability of the 1-twist state of the Langevin 
Eqs. ([I])-©. The open circle in (a) corresponds to the limit case D = C = studied in fit] ]. 

lengthy in this limit. 

We next compare the continuum Fokker-Planck description ([3]) with the discrete system, 
by means of Langevin simulations of Eqs. ([I])-©. To measure the stability of the 1-Twist 
state in Langevin simulations, we prepare the system in the 1-Twist state by randomly 
positioning the oscillators in space, and setting their phase to 9(x) = 2itx/L. We first let the 
system relax for 1000 units of time, so that the stationary shape of the angular distribution 
is reached. Then, we compute the ensemble average of the twist (m), where (...) denotes an 
average over 100 realizations of the numerical simulation, which we performed for different 
sets of parameter values C, D, and r/L. To evaluate the stability of the 1-Twist solution we 
set a threshold on the average twist: when (m) falls below 0.05, we take the 1-twist state to 
be unstable. Simulations were performed with iV = 1000 and L = 2tt. Other parameters are 
indicated in Fig. [5j The stability threshold in Fig. [5(a) was found by varying r for various 
fixed values of C, in Fig. [5(b) by varying D for fixed values of r, and in Fig. [5](c) by varying 
C for fixed values of D. We find a good agreement between the stability of the 1-twist 
state of the discrete system, and the stability of the 1-twist solution found by numerical 
integration of the Fokker-Planck equation. 
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VI. ATTRACTION BASINS 

Twisted solutions co-exist with global order and among themselves, Fig. HI As mobility 
is increased from low values, twisted solutions become unstable one by one, e.g. Fig. H(c) 
and Fig. [5(c), and global order is enhanced resulting in an increasing value of the ensemble 
average of the global order parameter as displayed in Fig. [2(b). Angular fluctuations can 
also destabilize twisted solutions, Fig. HJ However, while increasing angular fluctuations 
can decrease the number of twisted solutions competing with global order, the width of the 
angular distribution corresponding to global order also grows as a consequence of increasing 
angular fluctuations. The interplay between these two competing effects may be the cause for 
the non-monotonic behavior observed in the ensemble average of the global order parameter 
as angular fluctuations C increase, see green squares in Fig. [2(a). For these reasons, it is 
interesting to explore the size of the attraction basins of the different states that the discrete 
system exhibits. 

The fraction of realizations B(m) in which the system ends up in a particular state after 
a short time, starting from random initial conditions, is a measure of the size of the basins 
of attraction of the state. Coexistence of twisted states with up to m = 5 is shown in 
Fig. [6(a) for low mobility and phase fluctuations, corresponding to the bottom left corner 
of Fig. II^c) . As C is decreased, the attraction basin corresponding to global order shrinks, 
while the attraction basins of the m- twist states expand, Fig. [6(b). Decreasing C does 
not necessarily yield global order at the ensemble level because the global ordered state 
shares the phase space with the m-twist solutions: the observed value of (Z) results from 
the competition between an increase in local order for m = and a decrease of B(0), see 
Fig. [2(a). An increase in the mobility D leads to a contraction of the attraction basin of the 
m-twist states, in favor of global order, Fig. [6(c). When all m-twist solutions are unstable, 
5(0) = 1 and the global order state is the only attractor below C*, see also Fig. [2(b). 

VII. DISCUSSION 

We have investigated the effects of mobility in a generic ID model of locally coupled 
moving phase oscillators, and showed that oscillator mobility dramatically affects the collec- 
tive behavior of finite systems. More specifically, our results show that in low dimensional 
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FIG. 6: The size of attraction basins of m-twist solutions depends on angular fluctuations and 
mobility, (a) Probability B{m) that starting from a random initial condition the system ends up 
in a m-twist solution. B(0) and -B(l) as function of (b) C/7 and (c) D/j. Parameters as in Fig.0 

systems global synchronization is compromised by the presence of multiple m-twist states 
exhibiting local order. At the onset of local order, the system can fall into the global syn- 
chronized state. However, the coexistence of local order m-twist states implies that the 
attraction basin of the global synchronization state is reduced. Strong mobility of the oscil- 
lators destabilizes these m-twist states, and thus promotes global synchronization. 

In this paper we have considered a high density limit such that the connectivity of the 
system is never interrupted by gaps. In the dilute limit, gaps in the connectivity play a 
crucial role in the synchronization dynamics. This problem was studied in the context of 
moving neighborhood networks, and under the assumption of a fast exchange of neighbors, 
a mean-field condition for the existence and stability of the global synchronization state 
has been derived [23|. According to this study, whenever the global synchronized state is 



stable the system reaches global synchronization, regardless of its spatial dimensionality. In 
other words, the study overlooks the possibility of coexistence of multiple solutions. Our 
findings reveal a different role for mobility, unrelated to the existence and stability of the 
global synchronization state: mobility disrupts all these multiple solutions except for the 
global synchronized state. Extensions of the current study towards dilute systems will be 
the subject of further investigations. 



Two-dimensional systems display a similar phenomeno 



order states can now take other forms, e.g. vortexes 24J. It has been recently reported 



ogy, though the competing local 



that chaotic oscillators moving in a two-dimensional space can synchronize provided that 
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spatial dynamics is fast enough 25[. A related albeit different scenario occurs with chaotic 
advection mixing in two dimensional systems, where synchronization of excitable media is 
enhanced by strong mixing 

HQ. These results indicate that mobility may also enhance 
synchronization in two-dimensional systems. We speculate that global synchronization may 
be achieved by destabilizing local deffects, as we show here for ID systems. Further work is 
intended to clarify these issues. 

The theoretical framework introduced here may provide insight into other related prob- 
lems, as when movement is coupled to the oscillator phases. In this case synchronization 



can be interpreted as collective motion [28]. As a result of this coupling, strong spatial fluc- 
tuations and clustering effects dominate the system dynamics [29] , and global order prevails 
even in the thermodynamic limit 30]. 

Finally, a compelling biological application of our framework may be found in the verte- 
brate segmentation clock, where global coupling is a good effective description of the system 
because of the high mobility of cells [3jJ : by precluding the appearance of local defects, mo- 
bility promotes global synchronization. Moreover, it has been recently shown that mobility 
decreases the relaxation times to achieve synchronization in a model of the segmentation 



clock that allows for flipping between neighboring cells 
spatial patterns 



321 ] . However, this system also hosts 



and mobility is not accounted for in current distributed models. It 
will be interesting to see how mobility affects the synchrony recovery times and pattern 
reorganization after perturbation in such models. 
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